##SI Replication Code for Bush & Clayton "Facing Change" APSR 2022 

##set WD 

rm(list=ls(all=TRUE))
set.seed(1234)


library(lme4)
library(texreg)
library(effectsize)
library(stm)
library(quanteda)
library(tm)
library(geometry)
library(igraph)
library(Rtsne)
library(rsvd)
library(car)

setwd("...")

##APPENDIX B

lapop<-read.csv( "LAPOPData2.csv") #

#Table A2: Ideology as a moderating variable 

mixed2<-lmer(ClimStand ~ female + ideo + log(GDPPC) + (log(GDPPC) *female) + (1 +log(GDPPC) | country), data= lapop) #

mixed4<-lmer(ClimStand ~ female + ideo + log(GDPPC) + (log(GDPPC) *female) + (log(GDPPC) *female * ideo)  + (1 +log(GDPPC) | country), data= lapop) #
summary(mixed4) 

screenreg(list(mixed2, mixed4)) 

#Table A2, Pew data 

pew<-read.csv('PewData.csv')  #this dataset is in main directory 

mixedP2<-lmer(ClimStand ~ Female + Ideo + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data= pew) #
summary(mixedP2) #i 
screenreg(list(mixedP2))

mixedP4<-lmer(ClimStand ~ Female + Ideo + log(GDPPC) + (log(GDPPC) *Female) + (log(GDPPC) * Female * Ideo) + (1 +log(GDPPC) | Country), data= pew) #
summary(mixedP4) 

screenreg(list(mixedP2, mixedP4), digits=3)

##APPENDIX C

#Table A3: Negative gender gaps in low-income countries 

##only include countries w/ a "negative" gender gap 

modeldata<-read.csv("LAPOPRobust.csv")

neg<-subset(modeldata, country=="Guatemala" |country== "Bolivia" |country==  "Honduras" | country== "Haiti" | country== "Nicaragua" | country== "El Salvador" | country== "Peru" |country==  "Panama" | country== "Jamaica" |country==  "Ecuador" |country==  "St. Lucia" |country==  "Grenada")

lm1<-lm(ClimStand ~ female + country, data=neg) #
summary(lm1)

lm2<-lm(ClimStand ~ female + Edu + ideo + HHIncome + id_w_party + pol_interest + pol_know_rate + media + country, data=neg) #
summary(lm2)

texreg(list(lm1, lm2), digits=3)

##APPENDIX D

#Other surveys with climate questions 

#Figure A1: ISSP 

ISSP<-read.csv("DiffsISSP.csv")

cor.test(log(ISSP$GDPPC), ISSP $StandDiff) 

plot(ISSP $StandDiff ~log(ISSP $GDPPC), xlim=c(8,12), xlab="GDP per capita (logged)", ylab="Gender Gap (Female Mean - Male Mean)", main="Gender Gap in Perceived Seriousness of Climate Change \n by Level of Economic Development (ISSP Data)")
abline(h=0, col="black")
with(ISSP, text(StandDiff ~log(GDPPC), labels = X, pos = 4, cex=1))
text(6.8, 0.45, "R = 0.61, p < 0.001", cex=1, col="black")

#Figure A2: 2016 Pew

CountryLevel<-read.csv("Pew2016.csv")

plot(CountryLevel $Diffs ~log(CountryLevel $GDPPC), xlim=c(7,12), xlab="GDP per capita (logged)", ylab="Gender Gap (Female Mean - Male Mean)", main="Gender Gap in Perceived Threat of Climate Change \n by Level of Economic Development (Pew Data 2016)")
abline(h=0, col="red")
with(CountryLevel, text(Diffs ~log(GDPPC), labels = COUNTRY, pos = 4, cex=0.5))
text(7.6, 0.35, "R = 0.79, p < 0.001", cex=0.7, col="Blue")

#Figure A3: Arab Baromter  

Arab<-read.csv("ArabDiffs.csv")

Arab $GDPpc<-as.numeric(as.character(Arab $GDPpc))

plot(Arab $Diffs ~log(Arab $GDPpc), xlim=c(6,11), xlab="GDP per capita (logged)", ylab="Gender Gap (Female Mean - Male Mean)", main="Gender Gap in Perceived Seriousness of Climate Change \n by Level of Economic Development (Arab Barometer 2018 - 2019)")
abline(h=0, col="red")
with(Arab, text(Diffs ~log(GDPpc), labels = country, pos = 4, cex=0.7))
text(6.8, 0.3, "R = 0.69, p = 0.01", cex=0.8, col="Blue")

cor.test(Arab $GDPpc, Arab $Diffs) #corr = 0.68, p = =0.01

##APPENDIX E

##Excluding Anglo-phone and very poor countries (Tables A4 and A5)

#Table A4: AmericasBarometer 

modeldata<-lapop

NoAnglo<-subset(modeldata, modeldata$country != "Canada" & modeldata$country != "United States")
table(NoAnglo$country)

mixed1<-lmer(ClimStand ~ female + log(GDPPC) + (log(GDPPC) *female) + (1 +log(GDPPC) | country), data=
NoAnglo) #
summary(mixed1) #interaction term high significant (p<0.001)

mixed2<-lmer(ClimStand ~ female + ideo + log(GDPPC) + (log(GDPPC) *female) + (1 +log(GDPPC) | country), data=
NoAnglo) #
summary(mixed2)

mixed3<-lmer(ClimStand ~ female + HHIncome + Edu + ideo + log(GDPPC) + (log(GDPPC) *female) + (1 +log(GDPPC) | country), data=
NoAnglo) #
summary(mixed3)

screenreg(list(mixed1, mixed2, mixed3)) 

NoPoor<-subset(modeldata, log(modeldata$GDPPC) > (mean(log(modeldata$GDPPC)) - sd(log(modeldata$GDPPC))))
table(NoPoor$country)

mixed4<-lmer(ClimStand ~ female + log(GDPPC) + (log(GDPPC) *female) + (1 +log(GDPPC) | country), data=
NoPoor) #
summary(mixed4) #interaction term high significant (p<0.001)

mixed5<-lmer(ClimStand ~ female + ideo + log(GDPPC) + (log(GDPPC) *female) + (1 +log(GDPPC) | country), data=
NoPoor) #
summary(mixed5)

mixed6<-lmer(ClimStand ~ female + HHIncome + Edu + ideo + log(GDPPC) + (log(GDPPC) *female) + (1 +log(GDPPC) | country), data=
NoPoor) #
summary(mixed6)

screenreg(list(mixed4, mixed5, mixed6)) 

#Table A5: Pew  

modeldata2<-pew

NoAnglo<-subset(modeldata2, modeldata2$Country != "Australia" & modeldata2$Country != "United States" & modeldata2$Country != "Canada")
table(NoAnglo$Country)

NoPoor<-subset(modeldata2, log(modeldata2$GDPPC) > (mean(log(modeldata2$GDPPC)) - sd(log(modeldata2$GDPPC))))
table(NoPoor$Country)


mixedP1<-lmer(ClimStand ~ Female + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data= NoAnglo) #
summary(mixedP1) #interaction term high significant (p<0.001)

mixedP2<-lmer(ClimStand ~ Female + Ideo + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data= NoAnglo) #
summary(mixedP2) #including Personal Harm reduces coefficient by half. 
screenreg(list(mixedP2))

mixedP3<-lmer(ClimStand ~ Female + Edu + StandIncome + Ideo + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data= NoAnglo) #
summary(mixedP3) 

screenreg(list(mixedP1, mixedP2, mixedP3))

mixedP4<-lmer(ClimStand ~ Female + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data= NoPoor) #
summary(mixedP4) 

mixedP5<-lmer(ClimStand ~ Female + Ideo + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data= NoPoor) #
summary(mixedP5) #including Personal Harm reduces coefficient by half. 

mixedP6<-lmer(ClimStand ~ Female + Edu + StandIncome + Ideo + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data= NoPoor) #
summary(mixedP6) 

screenreg(list(mixedP1, mixedP2, mixedP3, mixedP4, mixedP5, mixedP6))

##APPENDIX F

##Additional robustness checks 
##Using Pew data

##Table A6

modeldata2 <- pew

#Merging in NEW contextual data
contextual.data<-read.csv('Pew_Contextual.csv')

combo.data<-merge(modeldata2, contextual.data, by.x="Country", by.y="country", all.x = TRUE, all.y = TRUE)

mixedP4<-lmer(ClimStand ~ Female + Edu + StandIncome +Tech + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data=modeldata2) #
summary(mixedP4) 

mixedP5<-lmer(ClimStand ~ Female + Edu + StandIncome +Tech + log(GDPPC) + (log(GDPPC) *Female) + climate_vulnerability2015 +
                (1 +log(GDPPC) + climate_vulnerability2015 | Country) , data=combo.data) #
summary(mixedP5) 

mixedP6<-lmer(ClimStand ~ Female + Edu + StandIncome +Tech + log(GDPPC) + (log(GDPPC) *Female) 
              + left_ideo + climate_vulnerability2015 + fossil_consump_2014 + hdi_2015 + 
                (1 +log(GDPPC) + left_ideo + climate_vulnerability2015 + fossil_consump_2014 + hdi_2015  | Country) , data=combo.data) #
summary(mixedP6)

##fully saturated model
mixedP7<-lmer(ClimStand ~ Female + Edu + StandIncome +Tech + log(GDPPC) + (log(GDPPC) *Female)  + left_ideo + climate_vulnerability2015 + oilrents_per_gdp + hdi_2015  + abs(latititude) + climate_HQs + (women_secondary * Female) + women_secondary + (1 +log(GDPPC) | Country) , data=combo.data) #
summary(mixedP7)

screenreg(list(mixedP4, mixedP5, mixedP6, mixedP7), digits=3)

##Table A7 (Americas Data) 

##fully saturated model w/ controls 

modelata<-read.csv('LAPOPRobust.csv')

mixed6<-lmer(ClimStand ~ female + HHIncome + Edu + ideo + pol_know_rate + id_w_party + pol_interest + married + n_child + risk_climate + media+ age + redist + urban+ log(GDPPC) + (log(GDPPC) *female) + (Edu*female) + (1 +log(GDPPC) | country), data=modeldata) #

screenreg(list(mixed6), digits=3)

##APPENDIX G

##Alternative dependent variable 
# Table A8 (Piew data)

mixed1<-lmer(SupportParis ~ Female + log(GDPPC) + (log(GDPPC) *Female) + (1 +log(GDPPC) | Country), data= combo.data) #
summary(mixed1) 

screenreg(list(mixed1))

##APPENDIX H

##Placebo issue: immigration  
#Figure A7

#left panel
plot(CountryLevel $Diffs ~log(CountryLevel $GDPPC), xlim=c(7,12), xlab="GDP per capita (logged)", ylab="Gender Gap (Female Mean - Male Mean)", main="Gender Gap in Perceived Threat of Climate Change \n by Level of Economic Development (Pew Data 2016)")
abline(h=0, col="red")
with(CountryLevel, text(Diffs ~log(GDPPC), labels = COUNTRY, pos = 4, cex=0.5))
text(7.6, 0.35, "R = 0.79, p < 0.001", cex=0.7, col="Blue")

#right panel
plot(CountryLevel $DiffsI ~log(CountryLevel $GDPPC), xlim=c(7,12), xlab="GDP per capita (logged)", ylab="Gender gap in immigration attitudes (standardized)", main="Gender Gap in Immigration Attitudes")
abline(h=0, col="red")
with(CountryLevel, text(DiffsI ~log(GDPPC), labels = COUNTRY, pos = 4, cex=0.5))
text(7.5, 0.13, "R = 0.12, p = 0.54", cex=0.7, col="Blue")


##APPENDIX M

##Netquest data from 2021 surveys in USA and Peru 

modeldata<-read.csv("USANetquest2021.csv")

##Table A9

model1<-lm(climateserious ~ meat + female + ptid1+ leftright_1 , data=modeldata) # meat consumption highly stat sig. 
summary(model1)

model2<-lm(climateserious ~ vehicle + female + ptid1 + leftright_1 + area, data= modeldata) #
summary(model2) 

model3<-lm(climateserious ~ drivetime + female + ptid1 + leftright_1 + area, data= modeldata) #
summary(model3) 

screenreg(list(model1, model2, model3), digits=3, stars = c(0.01, 0.05, 0.10))


#Table A10
#for men and women separately 

women<-subset(modeldata, female==1)
men<-subset(modeldata, female==0)

#for women and men separately

model1m<-lm(climateserious ~ meat+ ptid1+ leftright_1, data= men) #not 
summary(model1m)

model2m<-lm(climateserious ~ vehicle + ptid1 + leftright_1+ area, data= men) #
summary(model2m) 

model3m<-lm(climateserious ~ drivetime + ptid1 + leftright_1+ area, data= men) #
summary(model3m) 

model1w<-lm(climateserious ~ meat+ ptid1 + leftright_1, data= women) #not sig for women 
summary(model1w)

model2w<-lm(climateserious ~ vehicle + ptid1 + leftright_1 + area, data= women) #
summary(model2w) 

model3w<-lm(climateserious ~ drivetime + ptid1 + leftright_1+ area, data= women) #
summary(model3w) 

screenreg(list(model1m, model2m, model3m, model1w, model2w, model3w), digits=3, stars = c(0.01, 0.05, 0.10))

#Table A11

modeldata<-read.csv("PeruNetquest2021.csv")

model1<-lm(climateserious ~ meat + female + ptid + leftright_1, data=modeldata) # meat consumption highly stat sig. 
summary(model1)

model2<-lm(climateserious ~ vehicle + female + ptid + leftright_1 + area, data= modeldata) #
summary(model2) 

model3<-lm(climateserious ~ drivetime + female + ptid + leftright_1 + area, data= modeldata) #
summary(model3) 

screenreg(list(model1, model2, model3), digits=3, stars = c(0.01, 0.05, 0.10))


#Table A12 (American respondents)

modeldata<-read.csv("USANetquest2021.csv")

women<-subset(modeldata, female==1)
men<-subset(modeldata, female==0)

##SOCIOTROPIC COSTS  

modelM<-lm(climateserious ~ masc_fem_close_men_1, data= men) 
summary(model2)

modelW<-lm(climateserious ~ masc_fem_close_women_1, data= women) 

modelB<-lm(climateserious ~ close_gender + female + I(female * close_gender), data= modeldata) 
summary(modelB) 

screenreg(list(modelM, modelW, modelB), digits=3, stars = c(0.01, 0.05, 0.10))


#Table A13 (Peruvian respondents)

modeldata<-read.csv("PeruNetquest2021.csv")

women<-subset(modeldata, female==1)
men<-subset(modeldata, female==0)

modelM<-lm(climateserious ~ masc_fem_close_men_1, data= men) 
summary(model2)

modelW<-lm(climateserious ~ masc_fem_close_women_1, data= women) 
summary(model2)

modelB<-lm(climateserious ~ close_gender + female + I(female * close_gender), data= modeldata) #

screenreg(list(modelM, modelW, modelB), digits=3, stars = c(0.01, 0.05, 0.10))



